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Abstract 

We investigate features of Gravitational Waves (GWs) induced by primordial density fluctuations 
with a large amplitude peak associated with formation of Primordial Black Holes (PBHs). It is 
shown that the spectrum of induced GW is insensitive to the width of the peak in wavenumber 
space provided it is below a certain value, but the amplitude of the spectrum reduces at the peak 
frequency and decreases faster at low frequencies for a larger width. A correspondence between 
the GW amplitude and PBH abundance is also investigated incorporating the peak width. We 
find that PBHs with masses 10 20-26 g can be probed by space-based laser interferometers and 
atomic interferometers irrespective of whether the peak width is small or not. Further we obtain 
constraints on the abundance of the supermassive PBHs by comparing a low frequency tail of the 
GW spectrum with CMB observations. 
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I. INTRODUCTION 



Primordial black holes (PBHs) [l|-|3| are formed in the radiation dominated era if there 
exist density fluctuations of order of unity and overdensed regions enter the Hubble radius 
with masses determined by the horizon mass at that time 1 , 



47T 
~3~ 



ZPBH - M H = ^P^ 1 ) 3 , (1) 



together with those with smaller masses associated with the critical phenomena |i| whose 
contribution to the energy density is smaller. PBHs can therefore have a wide range of mass 
depending on the spectral shape of density fluctuations that generate them. 

If their mass is smaller than 10 15 g, they have evaporated away by now due to Hawking 
radiation jEJ . Abundance of PBHs in this mass range has been constrained by considering 
effects of emitted particles on the big-bang nucleosynthesis (BBN) j6j and the gamma-ray 
background 0] etc (see also 

PBHs heavier than 10 15 g, on the other hand, have not evaporated yet and can remain 
at present. These remnants could leave signatures in the present universe as astronom- 
ical objects. In the mass range M PB h ~ 1O 2 M — 10 5 M Q , they may be observed as 
intermediate- mass black holes (IMBHs) [loj . which are considered to be the observed ul- 
traluminous X-ray sources. They can also provide an astrophysical candidate for dark mat- 
ter (DM). Though PBHs with some range of mass are observationally excluded from being 
the dominant component of dark matter, there are no strong constraints in the mass range 
Mpbh ~ 10 20 g - 10 26 g (1O~ 13 M - 1O" 7 M ) 0, [HI. The methods to constrain in this mass 
range have been little known 2 . 

In order to account for these objects by PBHs, their present energy density is necessary 



to be n PBn h 2 ~ 1(T 5 - 1(T 2 for IMBHs [10] and n PB uh 2 = 0.1 for DM-PBHs [lj. To 
produce such a large number of PBHs, it is necessary to realize density fluctuations whose 
power spectrum has a high peak on the corresponding scales 3 . Some inflationary models are 
known to realize such a peaked spectrum (llj . 

At second order in perturbation theory, density fluctuations, or their associated scalar 
metric perturbations (scalar modes) inevitably generate gravitational waves (GWs, tensor 



modes) through scalar-tensor mode couplings |l5Hl7l]. Amplitudes of density fluctuations 
required for a substantial density of PBHs are so large that the second-order generation 
of GWs from them may well exceed the first-order counterpart (l8| . and amplitudes of the 



induced GWs could be large enough to be observed by GW observations [19] . 



In this paper, extending our previous Letter [19| , we further investigate features of induced 



GWs as a probe of PBH abundance. The spectrum of induced GWs carries information on 
that of density fluctuations. Therefore induced GWs provide information on density fluctu- 



ations [20j, and PBH abundance indirectly [19|, |21[. It has a peak at the peak scale of the 



1 



1. 
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We take 

2 In Ref. 13( , it is discussed that pulsar timing observations can constrain the abundance of PBHs with a 
mass -Mpbh ~ 10 25 g. 

One can also consider power law spectra with a blue tilt to realize a high amplitude on the small scales 
corresponding to the PBH mass considered here. These spectra, however, are not adequate for our 
purpose because they produce a larger number of PBHs on smaller scales and lead to a contradiction with 
constraints on the abundance of PBHs with smalllcr masses. 
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power spectrum of density fluctuations [19], which corresponds to the mass of the resultant 
PBHs. The corresponding frequency is estimated to be 10~ 9 Hz — 10 _8 Hz for IMBHs and 
10 _3 Hz — 1Hz for DM-PBHs. Fortunately, there exist observational means to constrain GW 



amplitudes in these frequency bands: constraints from pulsar timing observations |22l . |23 
for IMBHs and future observations by the space-based interferometers, LISA [3], DECIGO 
[25l | . BBO [2^], as well as the atomic gravitational wave interferometric sensors (AGIS) [2?J 
for DM-PBHs. 

In the previous Letter (l9| . we have shown the amplitude of the induced GW is large 
enough to be detected by these observations assuming the spectrum of density fluctuations 
is of delta-function type. In this paper, we extend the analysis to a spectrum with a finite 
width 4 to see how information on the spectrum of density fluctuations is imprinted on that 
of induced GWs. In the following sections, we show induced GW is useful to investigate 
PBH abundance even in this case. 

The organization of this paper is as follows. In section [Til we introduce PBH formation 
and give a relation between PBH abundance and amplitudes of scalar modes. In section IHH 
we estimate energy density of the GWs induced by scalar modes that lead to PBH formation 
taking into account a peak width. In section IIVI we give constraints on PBH abundance by 
using induced GWs. In section [V] we give a summary of this paper. 



II. PBH FORMATION 



In this section, we quickly introduce PBH formation with emphasis on a relation to 
density fluctuations, 5 = Sp/p, or their associated scalar metric perturbation \1>, which is 
defined by a trace part of a spatial metric. Here, the perturbation 6, \1/ are defined in the 
comoving gauge and in the longitudinal gauge, respectively. For later convenience, we use \l/ 
instead of 5 as a scalar degree of freedom. \l/ can be expressed in terms of 8 as ^ = (3/2)5 
at horizon crossing 28 . 



We use the Press-Schechter method 29 



to obtain an analytic formula that gives PBH 
abundance in terms of an amplitude of the scalar metric perturbation \F PBHs are formed 
in the radiation dominated era if the amplitude of the perturbation exceeds the threshold 
\l/ c at horizon crossing 5 . The fraction of the energy density of the Universe collapsing into 
PBHs with a mass M at the time of their formation, ppBH"(^f)/ptot, can be written in terms 
of the probability distribution of smoothed perturbation, Pr m as 6 



/3(M) 



Ppbh(M) 



Ptot 



"2 / Mn M 9 -^T- (2) 



dM 



Though the value of the threshold has not been determined definitely |3CH32| , we simply set 
^ c = 1/2 (S c = 1/3) {30} because our final result does not change much even if we employ 



21] 



4 See also Ref. 

5 In the followings, we assume all scales considered here re-enters the horizon in the radiation dominated 
era. This is satisfied if the reheating temperatura^isJiig her than 10 4 TeV 

6 The abundance /3(Mpbh) used in the literatures jalSl corresponds to an integral of our f3(M) with a lower 
limit, AipBH ■ 
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another value of \l/ c . The smoothed density fluctuation, ^r m , is defined by, 

# Rm {x) = J dV W(\x' - x\/R M )*(xf), (3) 

where W{x/Rm) is a window function, and the smoothing scale, Rm, is the scale of scalar 
modes that generate PBHs with a mass, M. From Eq.flTJ and clRm = H' 1 at horizon 
crossing, we can show Rm oc M 1 / 2 in the radiation dominated era, where H oc a~ 2 . Here, a 
is the scale factor of the Universe. 

We assume the probability distribution, Pr m , to be approximated by the Gaussian distri- 
bution because corrections from higher-order statistics of perturbations have been shown to 



be negligible in some inflationary models [34j, |35| . Then, Pr m is the Gaussian distribution 
with a variance, 

4 M = Jd\nkW(kR M ) 2 V*(k), (4) 

where W(kRM) is volume-normalized Fourier transform of the window function W(x/Rm), 
which is chosen to be a top-hat function here: W(JvRm) = 1 for kRu < 1 and WikRu) = 0, 
otherwise. V^(k) = (/c 3 /27r 2 ) (|\l/fc| 2 ) is a power spectrum of ^> at horizon crossing, which is 
assumed to have a peak at k p here. As \l/ is constant on super horizon scales (ggj, it coincides 
with a primordial one. 

Provided a power spectrum Vy, /3(M) can be estimated as Eq. (lA3p in Appendix IA~1 

/3(M)dM * V*{R£) —p-i exp f-^m d[ln(M/M p )], (5) 



' R M V K M ■ 

where we have introduced M p , which is the horizon mass when the peak scale, k~ l , enters 
the Hubble radius, k p = aH: 



M p = M eq 



k p \ i 5^ 



^eq/ \fl'*,eq 



-1/6 



K \ 2 / g*, P y 1 ? 6 



Here, M eq = 7 x 10 50 g, k eq = O.OlMpc -1 [14j, and #* >cq = 3.36 |36j is the horizon mass, the 



horizon-crossing wavenumber, and the effective number of relativistic degrees of freedom 
at the matter radiation equality, respectively, g*^ is the effective number of relativistic 
degrees of freedom at the time of PBH formation, which is ~ 10 for IMBHs (Mpbh ~ 
1O 2 M - 10 5 M G ) and ~ 10 2 for DM-PBHs (M PBH ~ 10 20 g - 10 26 g) 7 

As ppbh oc a~ 3 and PBHs are produced in the radiation dominated era, the density 
parameter of PBHs with a mass Mpbh at present, ^pbh(Mpbh), can be written in terms of 
/3(M PBH ) as 

/ fur \ -V2 / \ -1/3 
fipBH(MpBH)^ = /3(M PBH ) " tt m h 2 PBH / 9*, P 



M eq J V 9* 

-1/2 



The temperature at the horizon crossing of the peak scale is evaluated to be ~ 6McV — 60MeV for IMBHs 
and - 5TeV - 5 x 10 3 TcV for DM-PBHs, respectively. 
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where Q m h 2 = 0.13 [14| is the density parameter of matter. 

The factor in Eq.Q indicates that the PBH abundance is peaked at M p , which 
corresponds to the peak wavenumber of scalar modes, k p . Strictly speaking, however, the 
abundance has its maximum at Mpbh smaller than M p due to the exponential part in Eq.([5J). 
This is because the variance erjl is a monotonically decreasing function of M, so that the 
exponential part has a larger value for smaller M. Suppose that the power spectrum V^{k) 
has a flat peak extended over a wavenumber range k p e~ A < k < k p e A . If this region contain 
the most part of the total power of the spectrum, a\ =0 = A 2 , the variance cr^ , or the 
exponential part becomes constant for Rm < k p 1 e~ A , which corresponds to M < M p e~ 2A . 
In this case, the peak mass does not shift to a value smaller than M p e~ 2A , and is located at 
M p e~ 2A < Mpbh < M p . Therefore, a typical mass of PBHs corresponds to the peak scale of 
scalar modes, k~ l , within an accuracy of A. Power spectra of typical models such as those 



studied in [10|, |35[ sharply decrease for k > k p , and there is a one-to-one correspondence 
between the peak mass and the peak scale of scalar modes as Eq.([6]). Hereafter, we consider 
this type of spectrum. 

The total abundance of PBHs is given by an integral of Eq.(J7]) over M with an integration 
range 1O 2 M < M < 1O 5 M for IMBHs and M > 10 15 g for DM-PBHs. First, we consider 
the case that the peak of the power spectrum is steep enough so that the abundance has a 
peak at Mpbh — M p . When the peak mass, M p , is in the integration range, we can estimate 
the total abundance as 

W 2 = 4 x 10 13 % e (-^») ( J^-) (-*£-) " /3 , (8) 
A V 10 °g/ V 106.75/ ' v ; 

where we have used Vq(k p )2A ~ A 2 . 

In order to produce the PBHs that can explain the present IMBHs or DM, therefore, 
(\& c / ' v 4)e^~* c / 2 - 4 > must take a value between ~ 10 -14 and 10~ 8 . Hence, A 2 must take a 
value between 4 x 10~ 3 and 1 x 10~ 2 , which is larger than the observed value on CMB scales 
by a factor of O(10 7 ~ 8 ) [13]. In spite of a large difference in relevant values of the factor, 



M p 1/2 , in Eq.flSD, A 2 is of the same order of magnitude because it depends on this factor 
only logarithmically. 

The integration of Eq.([7j) is difficult to perform analytically for a general form of power 
spectrum. Here, we estimate the abundance for a top-hat type function spectrum, 

( A 2 

V ( J K) = < 2A f ° r l 111 ^/^' < A ' (9) 
otherwise, 

to understand how the expression (jSJ) of the total abundance is modified for large A. The 
abundance (j7j) is peaked at M PB h — M p e~ 2A due to the exponential part of /3(M PB h) and 
the factor M p ^. For large A with A > ^ 2 /4A 2 , the abundance ((7|) changes mainly by the 
factor Mpgy 2 at the peak. Hence, when the peak mass, M p e~ 2A , is in the integration range, 
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we can estimate the total abundance 8 as 



,13 *£,-*2/a*) / M p e^ 2A V 1/2 / x-i/a 



AA V !0 20 g / V 106.75 



(10) 



For large A, the upper limit of the integral in the first line should be replaced by 10 5 M Q 
for IMBHs. However, this replacement is irrelevant in the estimation because the integrand 
becomes exponentially small as mass increases. The total abundance (TTTJT) has the suppression 
factor (2A) _1 compared to Eq.(JH]). Hence, the required value of A 2 is larger than that for 
small A, though the increase is very slight because it depends on A only logarithmically. 
Hence, the power A 2 required to produce a relevant number of PBHs does not depend much 
on a detail shape of power spectrum in contrast to the peak mass. 

In the estimation above, we have not taken into account the critical phenomena, which 
leads to formation of PBHs with smaller masses corresponding to an amplitude of scalar 
modes as Mpbh = KMn(^> — ^c) 7 . If we employ this mass formula, PBH abundance has 



a peak at M PBH ~ 0.8M H - 1.1M H with a peak width, AM PB} JM U ~ 0(1) |4j, |37fl for 
A 2 ~ 4 x 10~ 3 - 1 x 10 -2 with K = 3.48 and 7 = 0.357 9,10 [Hi] when we use a delta 
function as V^. Therefore, there still exists the one-to-one correspondence between Mpbh 
and k p to within an order of magnitude if the critical phenomena are taken into account. 
This correspondence is also satisfied when a peak width is considered j33[. Further, the PBH 
abundance changes at most by a factor of O(10 _1 ). Hence, the corresponding values of A 2 
change very slightly thanks to its logarithmic dependence on the abundance. 

The scalar metric perturbations with large amplitudes as discussed above are expected 
to produce GWs with a much larger amplitude than usually expected through tensor-scalar 
mode couplings. In the following section, we estimate amplitudes of these induced GWs. 



III. GRAVITATIONAL WAVES INDUCED BY SCALAR MODES 
A. second-order gravitational waves 

Here, we give a brief review on generation of GWs from scalar modes as a second-order 



effect based on [16|, lljj [19|- We write a perturbed metric as 



ds 2 = a^) 2 [-e^drj 2 + e~™ g^dx 1 + V'drj)^ + V'drf)] (g tJ = 5 l3 + hy), (11) 
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9 
10 



In the case of IMBHs, correctly, we should use a total number density of PBHs instead of the energy 
density for estimating their abundance when their mass function has a width. However, the values of A 2 
required to explain IMBHs hardly change even if wc estimate it from the energy density (JTUJ) assuming a 
monochromatic mass function because PBHs with the peak mass mainly contribute to the total number 
density. 

Because we use \& instead of 5, the value of K referred here corresponds to if(2/3) 7 in [32| . 



32j | . the threshold is estimated to be ^ = 0.68 (S c = 0.45). Though we have used this value to estimate 



In 

the PBH abundance here, the abundance is almost the same even if we employ ^ c = 0.5. 
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where h\ = 0, and a gauge is chosen such that diV 1 = and <9j/i* = 0. Here, h l j = 5 lk hkj- 
The functions $ and ^ are the scalar modes and the functions V 1 represent vector modes. 
The tensor modes (gravitational waves) are represented by the functions h 1 -. ^ coincides 
with that used in the previous section at linear order. 

As we are considering second-order GWs, the tensor modes are second-order quantities 11 . 
Further, we assume the vector modes are also of second order because, at linear order, they 
only have a decaying mode and are exponentially suppressed in the standard inflationary 
cosmology. Expanding the Einstein equation up to second order, we obtain the evolution 
equation for induced GWs, 



ti" + ZHh)' - &% = (12) 

where a prime denotes differentiation with respect to the conformal time, rj, and % = a' /a. 
Here, the source term reads 

s r s = -2Vd r d s v + 4 dr(v + ?rV)d s (* + tt 1 *'), (13) 

3(1 "T Wj 

and V l r S j is a projection operator to a transverse, traceless part and w = p/p is the equation- 
of-state parameter. In practice, only the radiation dominated era with w — 1/3 is relevant, 
but we keep it as it is for future convenience. In the following analysis, we neglect anisotropic 



stress [17|, |38[, and set $ = \1/ at linear order. The anisotropic stress is expected to give only 
a small correction for GWs with a wavenumber k ~ k p because the peak scale, k' 1 , crosses 
the sound horizon before neutrino decoupling and the generation of the GWs mostly occurs 
at that time as explained later. In order to calculate the induced GWs up to second order, 
it is sufficient to use linear scalar modes. Hence, we only need to solve the linear evolution 



equation [28 



+ 6(1 + W) -V - wdH = 0, (14) 
l + 3w; 77 V ; 



for scalar modes. The evolution of the induced GWs is obtained by solving Eq.( fl2l with 

Then, working in the Fourier space, the solution of Eq.f fT2|) can be formally written as 

2 P 

hk(v) = J d?7 g k (r];r])a(r])VS k (fj), (15) 
where g k is defined as the solution for the following equation, 



a" 



9k+[k 2 --)gk = 6(r)-fj). (16) 

In other words, g k is the Green function of the evolution equation for ahk- It can be expressed 
by using a step function 9(rj — fj) as 5^(77; fj) = 9(r] — fj) sin[A;(?7 — fj)]/k during the radiation 
dominated era, where a oc 77. We define Fourier transform of tensor modes h k as 

j^m^- x M(vH(k) + hUvKM , (17) 



11 Such a treatment is justified when the induced GWs well exceed the first order counterpart. This condition 
is satisfied unless the first order GWs are amplified on the relevant scales. 
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where e^(fc), e*(fc) are the polarization tensors, which are normalized as 
£V . e ij(k)efj (— fc) = 2S rs . The function PS*, is defined in the same way. The pro- 
jection operator V projects k r k s in the Fourier transform of S r s to k 2 (l — ji 2 ) cos(20) for + 
mode and k 2 (l — fi 2 ) sin(20) for x mode where /i = k ■ k/kk and <fi is azimuthal angle of k r . 

The transfer function for Dk(r]), is defined by ^k(v) = -Dfc(^)^fc(0), which can be 
expressed by using a i^-th spherical Bessel function j„ as 

D k (rj) = (2v + \)\\x- v j v {x), x = ^k V , (18) 

where v = 1+ 2 3 ^ for an arbitrary value of w, and 



^) = -3 



x 2 



sin x 

cos(x) 



fc77/V3, (19) 



in the radiation dominated era with w = 1/3. Asymptotically, it behaves as 

D k {rj) oc h f ° r ^ X ' (20) 

| r/ -1 ' -1 cos (y/wkr/ — (v + l)7r/2) for y/wkf] ^> 1. 

Hence, ignoring oscillations, aVSk behaves as 



JV for a/^^p 7 ? < 1, 
I ?7 y for y/wk v r] > 1, 



(21) 



for the peaked power spectra of scalar modes. 

These behaviors show \aVS\~ \ has a maximum value at a time when the peak scale crosses 
the sound horizon in the radiation dominated era and the production of induced GWs mostly 
occurs at that time. After the sound- horizon crossing, hk evolves as a -1 , which is expected 
when the source term is irrelevant to the evolution of GWs. In the matter dominated era, 
the source term remains at a constant value and can be important again. However, it is not 
the case for the scales considered here because these scales re-enter the sound horizon at a 
time much earlier than the matter-radiation equality and the source term decays away in 
the radiation dominated era 12 . 

The power spectrum for tensor modes is defined as, 

2tt 2 

(h r k (v)h s M) = -^VnikMk + k^, (22) 
where the superscript r, s denote the polarization modes. 



12 We can estimate the scales where this statement is valid. Estimating k 2 h ~ S at the sound-horizon 
crossing time and requiring k 2 h > S at the present time, we obtain a condition MpBn/M cq < (1 + z eq ) -2 
for modes k ~ k p where we have assumed the peak scale k~ x crossed the sound horizon in the radiation 
dominated era. Here, z cq = 3176 Q is the redshift at the matter-radiation equality, so that the condition 
is written as Mpbh < 10 44 g ~ 10 10 M Q . If this condition is satisfied, the source term is also irrelevant for 
modes k < k p because these scales re-enter the sound horizon later than the sound- horizon crossing of the 
peak scale. 
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When the source term can be neglected, the energy density of tensor modes can be written 



in terms of their power spectrum for modes well inside the horizon [39 



1 (k^ 2 



£WM) = 3 ( ^ J Ph(k,ri). (23) 

Substituting the solution (Tl5"i) to Eq. (l2"3"j) . the energy density of induced GWs can be written 
in terms of the power spectrum of scalar modes as 

2 / k \ 2 P P 

^gw(M) = 3 ( J / d Vi dV2 a (Vi)a(V2)gk(v,Vi)9k(v,V2)Sk(ViiV2)- (24) 

Here, 

r°° _ r 1 k 3 k 3 ~ 
S k (vi,V2)= / dk d/i- ^(i-v 2 ) 2 f(k,\k-k\, Vl )f(k,\k-k\,r) 2 )v*(k)v*(\k-k\), 

Jo J -i |k — k| d 

(25) 

where f(ki,k2,rj) is a function written in terms of the transfer function for scalar modes, 
D k {rj) as follows, 

f(h,k 2 ,r]) = 2D kl (r))D k M + 3(1 A + w) [D kl ( V ) + U' x D' ki ( V )] [D k2 (rj) + H^D'^ (rj)} . (26) 

If we use k\ = k,k 2 = |k— k| instead of /c, /i as integration variables, Eq.f l24p can be rewritten 

as 

^gw(M) = ^ ( If dkidh VvMVvfaXl - I(k,k u k 2 , V ) 2 , 

(27) 

where 

/v 
drj a{r})g k {r)\ fj)f(ki, k 2 , fj), (28) 

and /ii2 = k • (k — k)//ci/c 2 . 



B. Gravitational waves induced by scalar modes with a peak 

We can calculate energy density of induced GWs by using Eq.( )24|) for a given power 
spectrum of primordial scalar pertuabtions, V^. Here, we estimate the energy density of 
induced GWs for a power spectrum with a peak. 

Here, we consider a power spectrum with a steep peak, or more correctly a spectrum 
whose most part of the total power is localized near a peak scale as assumed in the section 
llTl To simplify the analysis, we approximate this spectrum by a top-hat function <Q with a 
width A, 

{A 2 
— for \\n(k/k p )\ < A, ^ 
otherwise, 

where k p is the peak wavenumber and A 2 corresponds to the total power of the spectrum 
introduced in the section [Til As mentioned in the previous section, the peak mass of the 
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resultant PBHs is located at Mpbh = M v e for the spectrum (1291) . In the limit of small 
A, this spectrum can be further approximated by a delta function with respect to ln(k), 



(30) 



where the peak mass is given by M p . 



C. Delta-function type power spectrum 

Before giving a result for the spectrum, , we first investigate the simpler case, 

= Substituting the power spectrum fl30|) to Eq. fT24"l) . we obtain the energy density 

of induced GWs as follows, 



where 



2A 4 



aH 



k 



k 



- 2 

W) W^'M) 2 , (31) 



Id(k/k pj k p ri) — —I(k,kp,k p ,T)), 



(32) 



and 8(1 — k/2k p ) is a step function, which appears as a consequence of the momentum 
conservation law 13 . Here, we denote values at the time r\ v = k~ x with a subscript p. 

As we mentioned, the source term affects the evolution of GWs mainly at sound-horizon 
crossing of each k-mode. Hence, the function Id is nearly time-independent after the sound- 
horizon crossing. 

In the radiation dominated era, the function can be written as 



I d (k/k p , k p ri) 



d(k p fj) k p fjsm[k(r] -fj)]f(k p ,k p ,fj), 



(33) 



where Dk p (r]) in f(k p , k p ,rj) is given by Eq. ffT9l) . 

Here, we note that the function 1^ is not time-independent for the wavenumber k 
2k p /\/3. On subhorizon scales, the function a{rj) f (k p , k p , rj) / a p behaves as 



27 



2k p r] 



1 — cos 



2k p rj 

7T 



(34) 



Therefore, the phase of and that of / interfere constructively in the radiation dominated 
era and resonance amplification occurs for \k — 2k p / \/?>\ri CCi <ti 1, where rj eq is the conformal 
time at the matter-radiation equality [l6j]. Since the function / evolves as 77 -1 in the radiation 
dominated era, amplitudes of these modes grow logarithmically. 

We show the spectrum of the function l 2 d and the energy density of induced GWs in 
FigsJTJ In the right panel, we normalize the energy density of induced GWs by A 4 Q ra d{v) ■ 
By employing this normalization, we can extract its time dependence from ^gw because 
^rad(?7)Mad(?7p) = {9* J 9*) 1/3 {a p H p /aH ) 2 with fi rad (?7p) ~ I. 



13 In general, GWs with k < Nk p are generated at iVth-order. Hence, GWs with k > 2k p are also produced 
by higher-order effects though their amplitudes are much suppressed. 
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FIG. 1: The spectrum of the function /J given by Eq. ([33|) (left) and the energy density of induced 
GWs for the power spectrum (f50j) (right). In the left panel, the function /J estimated at k p i] = 
1.0 x 10 2 (solid line) and k p r] = 1.0 x 10 3 (broken line) are plotted. As mentioned in the text, 
ignoring oscillations, /J is time-independent. A time averaged one is also plotted (thick solid line). 
In the right panel, the plotted energy density is time averaged and normalized by the .4 4 f2 ra d. 



The energy density Oqw^ 2 has a peak at a frequency /gw = k p / (\/o7r) and scales as f 2 
for low frequencies. A typical amplitude of the induced GWs at the peak frequency is given 
by 

A sharp peak is observed in the spectrum of the function Ij at k ~ 2k p / a/3 in FigsJTJ This 
is due to the resonance amplification. 



D. Top-hat type power spectrum 



Next, we investigate a more general case, = V^ R \ To make a correspondence to 
Eq.()3TT) clearer, we introduce variables k c = (k\ + /c2)/2 and 8 = (k\ — k^jk. Substituting 
the power spectrum to Eq.f l7T|) . the energy density of induced GWs can be written as 

rt™\k/k p , V ) = (J^j JJ^ d(\nk c )d5 J{kJk,5)w(kJk,S,T])Q < §l(k/k c ,r]), (36) 



where 



and 



2\2 



w(k c /k,5,r)) = (1 - 5 ) 



J(k c /k,5) = ^- 
k c 



(JlY 


-1 


\2kJ 





I(k,k 1 ,k 2 ,ri) 
I(k, k c , k c , rf) 



(37) 



(38) 



which is a Jacobian. The function w(k c /k,S,rj) is nearly time-independent at sufficiently 
large r], because it depends on time only through I(k,ki,k2,T)). The integration domains, 
Di, D 2 , are given by 



D! = {|5| < l,k c >k/2}, 

D 2 = {| Hh/k^l < A, | 1h(V*p)I < A l- 



(39) 
(40) 
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Here, a restriction on the domain Di is a consequence of the momentum conservation law. 

The resultant energy density of the induced GWs, , is shown in FigsfSJ We have 
depicted the spectrum for A = 0.0, 1.0 x 10~ 3 , 1.0 x 10 _1 , 1.0 in the left panel and the right 
panel shows their values at k = k p normalized by that for the delta-function type spectrum, 

^GW' as a function of A at k p rj = 1.0 x 10 3 . 

The left panel shows the spectrum is almost the same for small A, but the spectrum with 
A = 1.0 differs from the others qualitatively. It behaves as oc k 3 for small modes 14 , has a 
plateau over a wavenumber range | ln(k/k p )\ < A, and a cut off appears at kjk p — 6 — 3. 
For k/k p < 2 sinh A, the integration range of 5 is determined by Di, and an additional factor 
k appears from the Jacobian (1551) . In fact, we can observe that an increasing rate of the 
spectrum changes at k/k p ~ 2sinh(0.1) ~ 0.2 for A = 1.0 x 10 _1 . The plateau reflects a 
behavior of the power spectrum of scalar modes, Vy. The main contribution of the integral 
(|36|) comes from a domain near the point (<5, k c /k) ~ (0, 2/v3), which is included by D2 for 
the wavenumber range \ \n(k/k p )\ < A, for large A. Hence, when A is large enough, the 

integral becomes irrelevant to a value of A, and VL^^\k/k p ~ 1,77) decreases as A~ 2 . An 
amplitude at the plateau is roughly estimated to be ~ coA _2 v4gw, where ^gw is a value 
calculated for the delta-function type power spectrum with an amplitude, V\&(2k / 'y/3) . Here, 
c is a numerical factor of O(0.1). The fc-dependence at the plateau, thus, reflects that of 

The right panel of Figs ]2] shows that ^ GW is almost constant for A < 0.1 and can be 
approximated by for A < 0.2. An increase observed at A ~ 0.3 is due to the resonant 
amplification. The resonant amplification occurs only for \k — (ki + k 2 ) / V3\r] cci <C 1 and 
contributions from these modes become smaller for larger A. Therefore we do not incorporate 
it to obtain a bound on PBH abundance in the next section. The energy density, furthermore, 
decreases as oc A~ 2 for large A as expected. These qualitative behaviors are expected to be 
the same even when we use other peaked power spectra, though threshold values of A may 
differ from the present values. 

To summarize this section, we have shown that the energy density of induced GWs (I3T)|) 
is independent of A and uniquely determines the total power of scalar modes, A, for small A 
with A < 0.2, whereas the energy density depends on A as oc A~ 2 for large A with A > 0.5. 
In the latter case, we need to estimate not only an amplitude of the induced GWs, but also 
a peak width of their spectrum to determine A. The peak width can be estimated if edges 
of the plateau are observed. 

IV. THE INDUCED GRAVITATIONAL WAVES AS A PROBE FOR THE PBH 
ABUNDANCE 

In the previous sections, we have seen how quantities of PBHs, (f2pBH^ 2 , A^pbh), and those 
of induced GWs, (Ogw^ 2 ,/gw)) are determined by those of scalar modes, (A,k p ,A). In 
particular, there is a one-to-one correspondence between (f2 PB H^ 2 , -^pbh) and (flcwh 2 , few) 
when a peak width of the spectrum, A, is sufficiently small. Hence we can use induced GWs 
to investigate PBH abundance indirectly. Here, we give constraints on PBH abundance from 
GW observations. 



The power spectrum should decay as Idln^/dlnfej > 1.0 at k < k p e A for this behavior to be manifest. 
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FIG. 2: The energy density of the induced GWs for the power spectrum (|29p for A = 0.0, 1.0 x 
10~ 3 , 1.0 x 10 , 1.0 (left) and its value at k = k p normalized by that for the delta-function type 
spectrum as a function of A (right). Both are estimated at k p r] = 1.0 x 10 3 . 



Mpbh and few can be written in terms of k p as Eq.flS]) and /gw = k p f (ySn), respectively. 
Hence we obtain a relation, 

/ GW = 0.03Hz ( %^ (^2_) " 1/12 . (41) 
JUW V 10 °g/ V 106.75/ V ; 

Thus, the mass range of the PBHs that can be an origin of IMBHs, Mpbh ~ 1O 2 M — 1O 5 M , 
and DM, M PBH ~ 10 20 g - 10 26 g, correspond to / GW ~ 10- n Hz - 10~ 9 Hz and / Gw ~ 
10 _5 Hz — 10 _2 Hz, respectively. In each frequency band, there are observations or obser- 



vational constraints: constraints from pulsar timing observations [22|, |23[ for IMBHs and 



observations in the near future by space-based interferometers, such as LISA (24); DECIGO 



25|, and BBO [26J, as well as AGIS [27J for DM-PBHs. 

The pulsar timing observations are sensitive to GWs with f > l/T where T is their 
data span. Hence, GWs with / ~ 10 _9 Hz — 10 _8 Hz are detectable by one-to-ten-years 
observations of the pulsar timing. Moreover, since pulsars are observed once every few 
weeks, detectable GW frequencies are limited to / < 10~ 7 Hz. Therefore, by using the 
pulsar timing observations, we can probe the abundance of PBHs with masses 10 _2 M Q < 
Mpbh < 10 2 M Q (T/35 yr) 2 . Since the GW spectrum extends up to / = \/3/gw! twenty- 
years observations could detect the induced GWs associated with IMBH-PBHs with masses 
Mpbh ~ 10 2 M Q . 

The space-based interferometers are sensitive to GWs with / ~ 10~ 5 Hz — 1Hz, which 
covers the enire mass range of the PBHs that are allowed to be DM, Mpbh ~ 10 20 g — 
10 26 g: LISA will have its best sensitivity fl G wh 2 ~ 10" 11 at / ~ 10~ 2 Hz (M PBH ~ 10 21 g), 
DECIGO/BBO and ultimate-DECIGO are planned to have sensitivities Qcwh 2 ~ 10~ 13 
and Q GW h 2 ~ 10~ 17 at / ~ lO-^z (M PBH ~ 10 19 g), respectively. 

In Figj3l we have plotted the energy density (|35j) corresponding to (f2pBH^ 2 , Mpbh) = 
(1O~ 5 ,1O 2 M ) and (Q PBU h 2 , M PBH ) = (10~\ 10 20 g) with the limit by the pulsar timing 
observations and the planned sensitivities of the space-based interferometers as well as those 
of AGIS and LIGO. We have also plotted envelope curves of the peak amplitude, A GW , 
which depends on /gw logarithmically with Qp^h 2 fixed at 10~ 5 and 10 _1 . The energy 
density fl36|) for large A is also depicted in Figj3j The spectra depicted correspond to 
(fipBH/* 2 ,MpBH,A) = (1O- 5 ,1O 2 M ,1.O) and {Q PB nh 2 , M PBH , A) = (10" 1 , 10 20 g, 1.0). We 
have also plotted envelope curves of the peak amplitude with Qp B nh 2 fixed at 10~ 5 and 10 _1 . 
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This figure shows that future observations planned by the Parkes Pulsar Timing Array 



(PPTA) project and the Square Kilometer Array (SKA) project |22J have sensitivities high 
enough to detect the induced GWs associated with IMBH-PBHs when A is small. The 
energy density becomes smaller for larger A. The amplitude of scalar modes, however, 
cannot be large on the CMB scales and A cannot exceed ~ 14. Therefore, the energy 
density of GWs cannot be smaller than O.IAqw/A 2 ~ 10~ 10 , if we take Aqw ~ 10~ 7 which 
is required to form appreciable amount of PBHs. Therefore, amplitudes of the induced GWs 
exceed the sensitivity of SKA even when A is large. To detect the induced GWs associated 
with IMBH-PBHs, we need to observe pulsars for a long period, T ~ 35 yr(MpBH/lO 2 M ) 2 . 
As already mentioned, however, since the GW spectrum extends up to / = \/3/gw, they 
could be detected in half the time. Moreover, Ref. |43| has proposed that a limit flc^h 2 ^ 
10~ 7 can be obtained for / ~ 10" 12 Hz - 10 -9 Hz (M PBH ~ 1O 2 M - 10 8 M o ) by measuring 
rotational parameters of pulsars. Hence, it might be possible to constrain IMBH-PBH with 
larger masses when A is sufficiently small. 

Further, it is clear from Figj3] that the planned sensitivities of the space-based interfer- 
ometers are high enough to detect the induced GWs associated with DM-PBHs irrespective 
of whether A is small or not. To avoid overproduction of evaporating PBHs 15 , A cannot be 
larger than ~ 10, where the expected energy density is estimated to be O.IAqw/ A 2 ~ 10 -11 . 
Even if we consider a broader peak extending to near CMB scales, A ~ 32, the energy den- 
sity is estimated to be 0.1Agw/A 2 ~ 10~ 12 , which can still be probed by using correlation 
analysis. Figure J3] further shows that we can observe edges of the plateau for A ~ 1 by these 
interferometers, and determine the peak width, A, from the observed shape of the spectrum. 
In this case, we can estimate the corresponding power of scalar modes by using the relation 
f^(Gw)^ 2 — O.IAqw/A 2 - The space-based interferometer can therefore determine whether 
PBHs can constitute DM in the Universe or not. 

Before ending this section, we mention other observational constraints on the induced 
GWs. The energy density of GWs can be constrained by considering their contributions to 



that of the universe at the BBN [44j. This constraints can be expressed as 

rf=°° 

/ d(log f)n GW h 2 < 0.2a ad /* 2 (iVc ff - 3.0), (42) 
Jf=o 

where N e s = 4.4±1.5 (68% CL) [3] is the number of relativistic degrees of freedom expressed 
in units of the effective number of neutrino species. Hence, the bounds from the BBN is 
conservatively given by Oqw^ 2 A < 1 x 10~ 5 . The energy denisty of induced GWs satisfies 
this constraints. 

CMB observations provide constraints f2 GW /i 2 < 10 _14 (//10 _16 Hz)~ 2 for / ~ 3 x 
10 _18 Hz — 1 x 10 _16 Hz by considering contributions of GWs to the temperature anisotropy 
For small A, the energy density of induced GWs scales as f 2 for / cq = 2 x 10~ 17 Hz < 
/ < few and as f° for / < / cq 16 , hence it is estimated to be ~ 10~ 17 - 10~ 21 for IMBHs 
and ~ 10 -31 — 10~ 37 for DM-PBHs at the frequency band relevant to CMB observations. 
Taking larger A leads to a smaller amplitude. Therefore, CMB observations provide only 



15 Note that the PBH abundance ((5]) is peaked at Mpbh smaller than M p . 

16 Provided the source term is irrelevant to the evolution of GWs, the transfer function of the induced GWs 
scales as f~ 2 in the radiation dominated era and as /~ 4 in the matter dominated era. Hence, we obtain 
^GW oc f° for / < / cq . 
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weak constraints on the abundance of the PBHs considered here. However, they constrain 
the abundance of the supermassive PBHs. Though the induced GWs are affected by the 
anisotropic stress for these supermassive PBHs, the change of their amplitude is less than 1% 



17(| . Hence, neglecting the anisotropic stress, the constraint on the power of scalar modes 
from CMB observations is given by 17 

^^(S" (43) 

for small A with A < 0.2 for the spectrum V\, . The mass range where the relevant 
constraints can be obtained is, thus, estimated to be Mpbh ^ 1O 9 M . For large A with 
A > 5.0 for the spectrum V^ R \ on the other hand, the energy density scales as f 3 for 
/ < f p e~ A and the amplitude at the peak scale is suppressed by a factor 0.1/A 2 . Hence, for 
large A, the constraint is given by 

A 2 < lO^Ae-l^- 1 ' ) ( ~ 3/4 . (44) 

Vio 9 M y v ; 



This constraint can be stronger than the constraint (1431) for large A thanks to the plateau 
extending to f p e~ A . Even though a dependence on the spectral shape of scalar modes 
exists, CMB observations can stringently constrain the abundance of the supermassive PBHs 
because the PBH abundance, (3, depends on A 2 exponentially. 



V. SUMMARY 

In this paper, we have investigated features of induced GWs as a probe of abundance of 



PBHs, extending the result of our previous Letter [19|] to a spectrum with a finite width. 
Amplitudes of density fluctuations that lead to PBH formation are so large that GWs in- 
duced by them through tensor-scalar mode couplings have a much larger amplitude than 
usually expected. The amplitude and the peak frequency of induced GWs are directly con- 
nected to the abundance and the mass of PBHs within an accuracy of the peak width. 
Fortunately, there is a correspondence between interesting mass scales and frequency bands 
which can be probed observationally: mass scale of PBHs that can be an origin of IMBHs 
and DM can be probed by pulsar timing observations and observations by space-based inter- 
ferometers, respectively. We have calculated the energy density of induced GWs expected 
when a relevant number of PBHs are produced. By comparing the result to the observa- 
tions, we have found that future long-term pulsar timing observations could detect /exclude 
IMBH-PBHs with masses Mpbh ~ 1O 2 M and measurements of rotation parameters of pul- 
sars more massive ones when the peak width is sufficiently small. We have further shown 
that the planned space-based interferometers, such as LISA, DECIGO, and BBO can deter- 
mine whether PBHs are DM or not, irrespective of the peak width. We have also discussed 



17 For PBHs with Mpbh — M eq ~ 10 17 M Q , the evolution of induced GWs at / > /gw is affected by 
the constant source term in the radiation dominated era, which has not decayed away until the matter- 
radiation equality. A mass range where the constraint is modified by the constant source term is estimated 
to be A/pbh — (0.3 — f.0)M eq by using the similar reasoning in footnote [T2] Hence the constraint (|43| 
can be used in most part of PBH mass range. 



15 



35 



/lg)IOgiQ(Mp BH 
30 25 20 



15 





-4 - 




-6 ■ 




-8 - 


o 




i-H 

Ol 


-10 > 


o 










-12 : 


5 






n 


-14 - 






-16 - 




-18 - 



pulsar 



advanced LIGO \ 
AGIS (ground) 




LISA\ 




DECIGO/BBD 



ultimate-DECIGO 



AGIS (space) 



IMBH 



DM 





-1 

-2 

-3 <? 

O 

-4 | 
o 

-5 
-6 
-7 



-10 



-6 -4 -2 

//Hz)ogi ( 



FIG. 3: Energy density of induced GWs (solid lines and broken lines) with a limit by pulsar timing 
observations and planned sensitivities of space-based interferometers (dotted lines). Thick solid 
lines indicate the energy density with the parameters (n PBR h 2 , M PB h) = (1O _5 ,1O 2 M ) (left), 
(10 -1 ,10 20 g) (right) for sufficiently small A. We have depicted those at the peak frequency as 
thick broken lines for fipBH^- 2 — 10 -5 (below), 10 _1 (above). Energy densities with A = 1.0 are 
also shown by thin lines. In estimating them, we have assumed that <7* iP is approximately 10 for 
the IMBH mass scales and 10 2 for the DM mass scales. Sensitivities of the space-based interferom- 
eters are depicted by using 40J with the instrumental parameters used in [4l|] for DECIGO/BBO. 
Sensitivities of the ground-based interferometer, LIGO [42|] are also plotted as a reference. The 



right vertical axis represents the corresponding power of scalar modes with small A and g 
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constraints on the abundance of the supermassive PBHs by CMB observations, and found 
stringent constraints. 
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Appendix A: Derivation of Eg. ([5]) 

Here, we give a derivation of Eq.([5]). The distribution function, Pr m , depends on M only 
through ct Rm . Hence, Eq.(j2]) can be rewritten as 

KM) = -2/^d* HM 

d * s „ MS&g*. (A!) 

In order not to overproduce PBHs, <j Rm /^ c < 1 should be satisfied. In this case, the 
variation of Pr m = (2ira RM )~ 1 ^ 2 exp(— ty 2 /2a R ) is mainly determined by its exponential 
part, and we can further rewrite Eq. flAlj) as 

where we omit a subscript Rm of \l/ for brevity. 
The integral in Eq. (lA2j) is estimated to be 



-Rm 

■*c^K M |f=f e —Pr m \v=H! c 



where we have used a RM /^ c < 1. Hence, the abundance © is roughly estimated to be 

/3(M)dM ~ V^Rg)—^- exp f-^p-) d[ln(M/M p )], (A3) 

where we have introduced M p , which is the horizon mass when the peak scale enters the 
Hubble radius. 
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